%% Simulate supply curve
close all; clc;
clear;

%% model parameters
T = 50; % number of age classes
age = 1:T; age = age(:);% vector of age classes

% Allometric biomass conversion 
w_a = 300*(1-exp(-age));

mature = ones(T,1); % dummy variable for maturity

% Natural morality 
m_a = zeros(T,1); % initialize
m_a(1:10)=[.49; .41; .32; .26; .22; .19; .17; .15; .13; .12];
m_a(11:20)=.11;
m_a(21:end)=.1;

% Beverton-Holt recruitment parameters 
alpha = 8.1; 
k = 10000;

%% catchability q(age)
q = (1 + exp(0*(2.4-age))).^(-1); % no gear selection 
q_a = (1 + exp(8*(2.4-age))).^(-1); % model gear selection

figure('Name','Catchability')
hold on
plot(age,q,'LineWidth',2,'DisplayName','Case I: \sigma = 0')
plot(age,q_a,':','LineWidth',3,'DisplayName','Case II: \sigma = 8')
xlim([1 10])
xlabel 'a, age'
legend('show','Location','southeast')

%% Simulate backward-bending supply
E = linspace(0,100,5e+4); E = E(:); % SS effort
c = 8000;

% ============== Single Peak Example ===================================%
H_v1 = zeros(length(E),1); % SS harvest
price_v1 = zeros(length(E),1); % SS price

N0 = zeros(length(E),1); % SS recruitment
Na = zeros(length(E),T); % age = 1:T

%-------------%-----------%-----------%-----------------%
for i = 1:length(E)
    e = E(i);
    z_a = e*q + m_a; % total mortality
    
    sum_z = zeros(T,1); % cumulative mortality from 1 to a
    for a = 1:T
       sum_z(a) = sum(z_a(age <= a));        
    end
    
    N0(i) = alpha*k - k./((mature.*w_a)'*exp(-sum_z));
    Na(i,:) = N0(i)*exp(-sum_z');
    
    H_v1(i) = max(0.1,Na(i,:)*((e*q./z_a).*(1-exp(-z_a)).*w_a));

end
%-------------%-----------%-----------%-----------------%
price_v1 = c*E./H_v1; % compute SS price
clear N0 Na

% ================= Multi-Peak Example==================================%
H_v2 = zeros(length(E),1); % SS harvest

N0 = zeros(length(E),1); % SS recruitment
Na = zeros(length(E),T); % age = 1:T
%-------------%-----------%-----------%-----------------%
for i = 1:length(E)
    e = E(i);
    z_a = e*q_a + m_a; % total mortality
    
    sum_z = zeros(T,1); % cumulative mortality from 1 to a
    for a = 1:T
       sum_z(a) = sum(z_a(age <= a));        
    end
    
    N0(i) = alpha*k - k./((mature.*w_a)'*exp(-sum_z));
    Na(i,:) = N0(i)*exp(-sum_z');
    
    H_v2(i) = Na(i,:)*((e*q_a./z_a).*(1-exp(-z_a)).*w_a);

end
%-------------%-----------%-----------%-----------------%
price_v2 = c*E./H_v2; % compute SS price

%% Na/N0 of first and second peak
[Hmax1, idx1] = max(H_v2);
H_copy = H_v2(E>=4);
[Hmax2, idx2] = max(H_copy);

ratio_peak1 = Na(idx1,:)./N0(idx1);
ratio_peak2 = Na(idx2,:)./N0(idx2);

figure('Name','NaN0 ratio')
hold on
plot(age',ratio_peak1,'bo-','LineWidth',2,'DisplayName','first peak')
plot(age',ratio_peak2,'*:','LineWidth',3,'DisplayName','second peak')
xlim([1 10])
xlabel 'a, age'
ylabel 'ratio of N_a and N_0 at harvest peaks'
title 'Steady-state population by age class'
legend('show')
%% Plot Backward bending Supply
figure('Name','Supply')

subplot(2,1,1)
hold on
title '(a) two cases of backward-bending supplies'
plot((1e-3)*H_v1,sqrt(price_v1),'LineWidth',2, ...
    'DisplayName','Case I')
plot((1e-3)*H_v2,sqrt(price_v2),':','LineWidth',3, ...
    'DisplayName','Case II')
xlabel 'Catch Quantity (MT)'
ylabel 'Price (square root scale)'
ylim([0 1])
legend('show')

subplot(2,1,2)
hold on
title '(b) single-peak vs. multi-peak supplies with respect to effort'
plot(sqrt(E),(1e-3)*H_v1,'LineWidth',2, ...
    'DisplayName','Case I')
plot(sqrt(E),(1e-3)*H_v2,':','LineWidth',3, ...
    'DisplayName','Case II')
xlabel 'Steady-State Effort (square root scale)'
ylabel 'Catch Quantity (MT))'
legend('show')
